Happy Pi Day everyone!!

Since the date today approximates pi (3.14), I thought I would celebrate by deriving an approximation of pi using light sabors.

We start in a futuristic world where urban wildlife have aquired light sabors and taken over Montreal.

Future squirels battle for food scraps near homeless man's dome tent.

Future squirels battle for food scraps near homeless man’s dome tent.

To battle these savage beasts, we must learn to make our own light sabor.

set.seed(1.1)
light_sabor_length <- 20
number_of_animal_sabors_in_a_park <- 5000
distance_between_transects <- 40.1

point_1 <- c(-(light_sabor_length/2),(light_sabor_length/2))
point_2 <- c(0,0)
single_sabor <-cbind(point_1, point_2)
rownames(single_sabor)<- c("point_1","point_2")
colnames(single_sabor)<- c("x","y")
single_sabor
##           x y
## point_1 -10 0
## point_2  10 0
Represent our lightsabor as a line with no width. A single beam of light.

Represent our lightsabor as a line with no width. A single beam of light.

Then distribute the army of sabor-wielding critters across the landscape.

First we rotate around the origin to set the instantanious angle of each animal’s individual sabor.

  angle <- runif(1, min=1, max=360)*pi/180
  xy <- as.matrix(single_sabor)
  
  # Rotation
  cos.angle <- cos(angle)
  sin.angle <- sin(angle)
  after_rotation <- xy %*% t(matrix(c(cos.angle, sin.angle, -sin.angle, 
        cos.angle), 2, 2))
  
after_rotation
##              [,1]      [,2]
## point_1  1.100398 -9.939272
## point_2 -1.100398  9.939272
We rotate the lightsaber around the origin to change its angle.

We rotate the lightsaber around the origin to change its angle.

Then we move (a.k.a. translate) the sabor to a location on the landscape where the critter is located. We assume that all the critters are battling vigously and indescriminintly with each other so they are randomly distributed across space and the angles of all the sabors are uniformly distributed.

  # Translation
  translation.x <- runif(1, min=10, max=90)
  translation.y <- runif(1, min=10, max=90)
  translation <- matrix(
    c(translation.x,translation.x,translation.y,translation.y), 2, 2)
  new_coord <- after_rotation + translation
  colnames(new_coord) <- c("x","y")
  rownames(new_coord) <- c("Point_1","Point_2")
After the angle is set, translate the lightsaber to a random location on the battlefield.

After the angle is set, translate the lightsaber to a random location on the battlefield.

Package that process into a function:

critter_sabor <- function(){
  single_sabor <-cbind(c(-5,5),c(0,0))
  angle <- runif(1, min=1, max=360)*pi/180
  xy <- as.matrix(single_sabor)
  
  # Rotation
  cos.angle <- cos(angle)
  sin.angle <- sin(angle)
  after_rotation <- xy %*% matrix(c(cos.angle, sin.angle, -sin.angle, 
        cos.angle), 2, 2, byrow=TRUE )
  
  # Translation
  translation.x <- runif(1, min=10, max=90)
  translation.y <- runif(1, min=10, max=90)
  translation <- matrix(
    c(translation.x,translation.x,translation.y,translation.y), 2, 2)
  new_coord <- after_rotation + translation
  colnames(new_coord) <- c("x","y")
  rownames(new_coord) <- c("Point_1","Point_2")
  
  return(new_coord)   
}

critter_sabor()
##                x        y
## Point_1 21.93899 84.59096
## Point_2 30.33012 79.15139

Repeat it many times to create and army

for(i in 1:number_of_animal_sabors_in_a_park){
critter_sabor()
}
An army of lightsaber-wielding critters in the middle of La Fountaine park. This is the position of each critter's lightsaber at the moment we pass them on our transect.

An army of lightsaber-wielding critters in the middle of La Fountaine park. This is the position of each critter’s lightsaber at the moment we pass them on our transect.

surface_line <- function(line_height = 0){
  matrix(c(0,100,line_height,line_height), 2, 2, byrow=FALSE)
  }
##      [,1] [,2]
## [1,]    0   20
## [2,]  100   20
My transect through the army of critters.

My transect through the army of critters.

one_sabor <- critter_sabor()
one_surface_line <- surface_line(30)
line.line.intersection(one_sabor[1,], one_sabor[2,], one_surface_line[1,], 
                       one_surface_line[2,], interior.only = TRUE)
## [1] 75.73794 30.00000
Build a mechanism to identify when a lightsaber has crossed my transect line.

Build a mechanism to identify when a lightsaber has crossed my transect line.

Make a list of surface lines and a list of matches

list_of_matches <- array(NA,dim = c(2,2,number_of_animal_sabors_in_a_park))
list_of_surface_lines <- array(NA,dim = c(2,2,6))

for(i in 1:number_of_animal_sabors_in_a_park){
  list_of_matches[,,i] <- critter_sabor()
}

line_placement <- seq(0,100, by=20)
for(i in 1:6){
  list_of_surface_lines[,,i] <- surface_line(line_placement[i])
}
Me and the team preparing for our transects.

Me and the team preparing for our transects.

After the angle is set, translate the lightsaber to a random location on the battlefield.

After the angle is set, translate the lightsaber to a random location on the battlefield.

crosses <- matrix(NA, 1, 10, byrow=TRUE)

for(i in 1:dim(list_of_matches)[3]){
for(j in 1:dim(list_of_surface_lines)[3]){
  crosses <- rbind(crosses, c(
    line.line.intersection(list_of_matches[1,,i], list_of_matches[2,,i], 
                           list_of_surface_lines[1,,j],     
    list_of_surface_lines[2,,j], interior.only = TRUE), 
    
    list_of_matches[1,,i], list_of_matches[2,,i], list_of_surface_lines[1,,j], 
    list_of_surface_lines[2,,j])
    )
  }
}

crosses <- crosses[which(is.na(crosses[,1]) == FALSE),]
colnames(crosses) <- c("x.cross", "y.cross", "x.p1", "y.p1", "x.p2", 
                    "y.p2", "x.surface.p1", "y.surface.p1", "x.surface.p2", 
                    "y.surface.p2")
head(crosses)
##       x.cross y.cross     x.p1     y.p1     x.p2     y.p2 x.surface.p1
## [1,] 63.08661      20 64.26028 15.53502 61.71804 25.20647            0
## [2,] 14.37017      80 12.61642 83.36074 17.24275 74.49524            0
## [3,] 49.44207      80 44.92515 85.23335 51.45904 77.66313            0
## [4,] 47.22221      20 49.37535 17.30391 43.13500 25.11787            0
## [5,] 24.43356      20 23.31048 15.67796 25.82545 25.35654            0
## [6,] 56.12789      20 53.77911 24.62028 58.31078 15.70602            0
##      y.surface.p1 x.surface.p2 y.surface.p2
## [1,]           20          100           20
## [2,]           80          100           80
## [3,]           80          100           80
## [4,]           20          100           20
## [5,]           20          100           20
## [6,]           20          100           20
All the lightsabers we struck on my transect.

All the lightsabers we struck on my transect.

pi_estimate <- (2* light_sabor_length * number_of_animal_sabors_in_a_park)/(distance_between_transects * length(crosses[,1]))
options(digits=20)
pi_estimate
## [1] 3.1427417593382642735

We take the ratio of the lightsabers we made contact with againt the ones we missed and we get a great approximation of Pi!

Happy Pi Day!